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I. INTRODUCTION 

Frequency-shift-position (fsp) modulation sends a sample x e (0, T) 
by changing the frequency of a wave from coo to toi at the time x. We 
consider the problem of optimal mean-square estimation of x from an 
fsp signal to which white noise has been added in transmission. The 
best estimator, given by a known formula of nonlinear filtering, con- 
tains stochastic integrals and is hard to implement. An approximation, 
obtained by omitting 0(wf I ) terms, is readily implemented over the 
next interval (0, T) by ordinary differential equations driven by the 
observed signal. 

We pose the optimal least-squares demodulation problem for the 
following simple communications system: The signal to be sent con- 
sists of successive i.i.d. random variables x having a density b with 
support (0, T); the channel provides an fm wave, of random initial 
phase 9, which shifts frequency from <o to wi at the point x during 
(0, T); the resulting signal is 

_ (R COS(0 + OJot, t<x 

s '~ [# cos(0 + coo* + (t-x)ui), T>t>x; 

this signal is observed through white noise, the received signal or 
observation being y, given by 

dy, = St dt + db t , 

where b t is a Brownian motion independent of x; the process is 
repeated over each interval of length T, and the problem is to guess or 
calculate a good estimate of x from the observation {y,, < t < T). 
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This problem is resolved by using nonlinear filtering to give an exact 
but not recursive least-squares estimator x based on {y t , < t < T). 
The estimator x is then analyzed for ways of calculating it in real time. 
First, it is re- expressed in terms of ordinary integrals by eliminating 
the stochastic differentials. Second, we find that by ignoring two terms 
that are 0(i? 2 /min{a>o, coi}) we can use a Bessel integral and some 
simple differential equations to calculate an approximate version of x 
over the next interval of length T. This version will be very close to x 
when the carrier frequencies coo and on are large compared to the 
amplitude R. To make the integrals appearing in the likelihoods we 
will use suitably dimensionless, we let R have unit sec~ 1/2 . 



II. FORMULA FOR x 

At time T, the best mean-square estimate of x given the observation 
{y„ < s < T) is given by the Kallianpur-Striebel or Bayes formula 1 



H 



x = D~ x — d6 tf{t) dt exp\R \ cos(0 + w s) dy s 

R 2 f f T 

— — I cos 2 {8 + ioos) ds + R I cos(0 + oiot + coi(a — t)) dy u 

2 Jo Jt 

R 2 f T , 1 

— — I cos 2 (0 + uot + ui(u — t)) du>, 

where D = normalizer = same integral with f{t) for tf(t). We first 
reduce this formula to ordinary integrals. Since the signal s(t) is of 
bounded variation, the mixed quadratic variation ( ys ) is 0, and the 
differentiation formula justifies an integration by parts that expresses 
the exponent above as 

f' R 2 C 

Ry,cos(6 + u t) + Ruo j s sin(0 + co s) ds — — I cos 2 (6 + u s) ds 

+ Ry T cos(d + u t + «i(T - t) - Ry ,cos(d + toot) 

r 

+ Rcoi ,y u sin(0 + u t + ui(u - t)) du 

R 2 C 
— — cos 2 (0 + w £ + ui(u - t)) du. 

The stochastic integrals are thereby eliminated. Also, 
R cos(0 + aot + ui(T - t)) = s(T), 
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and the cos 2 integrals can be evaluated thus: 

6 + UQt 

cos 2 (0 + toos) ds = wo 1 | cos 2 u du 





{ 



-,1 , . , e+Unt t sin(20 + 2co t) - sin 26 

- «0 q \U + Sill U COS w)fl °' = - + - 

Z 2 4uo 

and similarly 

I cos 2 (0 + w * + wi(w - t)) du = tor 1 - (u + sin w cos u) e g i^ t l + ^ {T - l) 



2 



_ T - t | sin(2fl + 2coq< + 2coi(r - t)) - sin(2fl + 2w t) 
2 4toi 



We thus arrive at the formula for x: 

•T 



* = ° l 2^\ de \ W) dt eXp | R y TCOS W + <^ot+u l (T- t)) 
4- Rioo ;y s sin(0 + coos) ds + Run y s sin(d + u t 



, . ,,, , R 2 T n2 sin(20 + 2cc t) - sin 20 

+ oji(u - t)) du R* - 

4 8co 

_ 2 sin(2fl + 2u3 t + 2un(T - t)) - sin(2fl + 2ui t) \ 

8«i J' 

again with D the same with f(t) in place of tf(t). Note that the i? 2 T 
term can be ignored, since it occurs also in D. 

The d6 dt integration can be interchanged by absolute convergence, 
so we look at the dependence of the exponent on 6. This is of the form 



a, cos 6 + fitsin + 



\a>o A u>\J 



where with 

•T ~T 

A= J y s sinu)isds, B= I y s cos cois ds, co = <o — coi 



f • f 

I > s sin cois ds, B = 

Jo Jo 

= /?U5o I .VsSU 

Jo 



a, = Ruo y s sin u s ds + y T R cos(u t + coi(T — t)) 
Jo 

+ Rui(B sin u)t + A cos uot) 

— Ru\cos ut y ssin u\s ds — Rooisin at J y s cosuiisds 



'coisin u)t J 
Jo 
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fit = Rioo y s cos coos ds — yrR sin(to £ + wi(T — t)) 



+ R(ji\{B cos ut — A sin ut) 

— Rioicos ut y s cos wis ds + i?uisin ut I y s smuisds. 



y s sh 

Jo 

The vector (a, /?) satisfies the differential equations: 

— coiJtoo I Vs ( I wos ds. 

Jo \- sm / 

It can be expressed in terms of the quantities 



s{t) = ott- Ruo j'ssin 100S ds, 
Jo 

c(t) = fit - Roio I ^cos coos ds, 
Jo 

which satisfy the nicer system 

*\ f ° A( s \ d /sinwoi 
c) = »[-l o){c)- R " iy \co S » t 

The initial conditions for a, p, and also for s, c, are 
so = oto = RwiA + Ryrcos <oi T 
Co = po = RuiB — i?yrsin <oi T. 

III. APPROXIMATION FOR SMALL /?7(<*>o A u, ) 

The "double-phase" 0(R 2 /uo A wi ) term in the exponent is indepen- 
dent of y, and considerable simplification results if we ignore it. The 
numerator of x is then nearly 

I tf(t)— e a ' cose+li,sin6 dd dt, 

Jo 27T J-n 

and we can use the Bessel function integral representation 

— ( e'™' + " h " dd = I (Va 2 + jB 2 ) 

277 I 
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to "do" the dd integral approximately, to get 

isr 1 [ tf(t)I (Ja? + tf) dt, (1) 

Jo 

again with D the same integral with the t factor left out. 

IV. IMPLEMENTATION 

We can now use the differential equations for a t and fit to suggest 
how to realize the approximate estimator (1) in real time. Since the 
initial values ao and fi depend on both 

r 

A = y s sin u\s ds 



■I 



B = I ^ s cos wis ds 

Jo 

as well as on yr, we shall first wait till the end of the current interval 
of length T, all the while "beating" the received signal y s with the 
"second" fm frequency coi and integrating, to end up with A and B at 
time T. Noting the value of yr, we obtain a and fi and start to solve 
the differential equations for a t and fit during the next interval of 
length T. This is achieved by having, in the preceding interval, beaten 
y s with the "first" fm frequency coo, integrated, and passed the results 
through a delay of size T to obtain, at time T + t, each of 

o y '[-^) ao8d8 ' 

this vector being precisely what we need to add to the solution of the 
differential equations 

l)" a (-i o)(" )"***(-*)*' 

to get a, and fit. We also delay y, by T so as to have it available in real 
time as soon as So and Co are known, so as to use it to drive the system 
for s, c in the last equation. 

Once a t , fit are available, we can calculate io(Va? + fif), and thence 
the integrals making up the ratio (1) that approximates x. Figure 1 is 
a block diagram for this implementation. 

Since the appro ximate e stimator depends on y, only through the 
"magnitude" r = Va 2 + fi 2 , some simplifications occur upon a change 
to polar coordinates. Putting also <$> = arctan a/fi, we obtain these 
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.1 



J, 



1 294 THE BELL SYSTEM TECHNICAL JOURNAL, SEPTEMBER 1 980 



equivalent equations in terms of the new variables r = r, and <f> = fy: 



r = —uRytCos{<}> — uot) — 00R000 I _y s sin(co — ioqs) ds 



(f> = —icRy t sin(a) — coot) — ooRcoo I y s co&(4> ~ w os) ds. 

'o 

It is perhaps not surprising that these equations should bear some 
resemblance to those arising in the study of the pendulum, the phase- 
locked loop, and the FM-with-feedback circuit. 



V. BASEBAND VERSION 

Instead of shifting frequency from wo to wi at the sample value x e 
[0, T], we can imagine changing a signal from the constant value ro to 
ri at the sample value (time) x. This leads to a "baseband" version of 
the problem just studied. The signal is now 



St = 



To < t < X 

n x < t < T, 



and the observation is dy, = s, dt + db t again. The optimal mean- 
square estimate x is 

xD~ l J tf(t)e**-wi*+ntyr-»)-v*HT-t) dt 

Jo 

where the normalizer D is the same integral with f(t) for tf{t). Figure 
2 is an exact implementation of this estimator in real time. The terms 
in yr and r\T cancel with the same ones in D, so can be ignored. 
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Fig. 2— Baseband version of estimator. 
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